Treatment order is not a significant predictor of wingbeat frequency in the model with percent load (though it is very close, p ~ 0.07)
Percent load and treatment suggest very similar models in terms of wingbeat frequency
Treatment order is not a strong predictor of stroke amplitude in the model with percent load
wingbeat freq = size Met rate = mass carrying amplitude = percent loading
Susie’s (hypotheses): lm( frequency ~ ITspan + trt order)
lm( metabolism ~ total mass)
lm( stroke amplitude ~ %load + trt) (those are probably the same thing)
How heavy the bee is loaded overall, affects how much the amplitude and frequency change
Questions to look into: ** A really heavy loaded bee will have a smaller change in metabolism vs. light bee ** Suggests a non-linear relationship or random slopes that depend on bee size — the cost for additional weight goes down as you get heavier loaded. — not the case for the change in amplitude (change in load vs. change in amplitude is linear) — not the case for metabolism or wingbeat frequency.
** Hypotheses: 1. Wingbeat freq moves around for a number of reasons (tired, cold, added mass) 2. Stroke amplitude rescues wingbeat frequency 3. As you get really close to your max, adding more weight doesn’t seem to cost as much.
Analyses of respirometry data:
Install required packages and read in data Define custom function for evaluating VIF with multilevel models
ipak <- function(pkg){
new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
if(length(new.pkg)) install.packages(new.pkg, dependencies = TRUE)
sapply(pkg, require, character.only = TRUE)
}
packages <- c('ggplot2', 'car', 'lme4', 'gsheet', "MASS", 'influence.ME', 'sjPlot')
ipak(packages)
## ggplot2 car lme4 gsheet MASS
## TRUE TRUE TRUE TRUE TRUE
## influence.ME sjPlot
## TRUE TRUE
theme_set(theme_bw())
# read in data -- google sheet called "Bumble mumble grumble"
bdta <- gsheet2tbl("https://docs.google.com/spreadsheets/d/1GUUoFq41Ep3sJNFXiyMcp7fZhYmQCzTcLnhVXdr4WRo/edit?usp=sharing")
(used later)
vif.mer <- function (fit) {
## adapted from rms::vif
v <- vcov(fit)
nam <- names(fixef(fit))
## exclude intercepts
ns <- sum(1 * (nam == "Intercept" | nam == "(Intercept)"))
if (ns > 0) {
v <- v[-(1:ns), -(1:ns), drop = FALSE]
nam <- nam[-(1:ns)]
}
d <- diag(v)^0.5
v <- diag(solve(v/(d %o% d)))
names(v) <- nam
v
}
# create centered and squared variables
bdta <- within(bdta, {
# new variables
percLoad <- (Mstarved + load) / Mstarved
totalMass <- Mstarved + load})
bdta <- within(bdta, {
# centered variables
percLoad_cent <- scale(bdta$percLoad, center = TRUE, scale = FALSE)[,1]
totalMass_cent <- scale(totalMass, center = TRUE, scale = FALSE)[,1]
wbf_cent <- scale(wbf.aud., center = TRUE, scale = FALSE)[,1]
strokeAmp_cent <- scale(stroke.amplitude, center = TRUE, scale = FALSE)[,1]
wingVel_cent <- scale(wing.velocity, center = TRUE, scale = FALSE)[,1]
#squared
percLoad2 <- percLoad_cent^2
totalMass2 <- totalMass^2
wbf2 <- wbf_cent^2
strokeAmp2 <- strokeAmp_cent^2
wingVel2 <- wingVel_cent^2
# change factor properties
# convert trt order to factor
Treatment.order <- as.factor(as.character(Treatment.order))
# change reference level to unloaded
Treatment <- factor(Treatment, levels = c("UL", "L"))
})
summary(bdta)
## Bee.ID Treatment.order Treatment Mstarved
## Length:60 1:30 UL:30 Min. :0.0767
## Class :character 2:30 L :30 1st Qu.:0.1213
## Mode :character Median :0.1378
## Mean :0.1411
## 3rd Qu.:0.1670
## Max. :0.1909
## load IT.Span single.wing.area av.resp..CO2.mL.hr.
## Min. :0.01650 Min. :3.860 Min. :20.03 Min. : 5.507
## 1st Qu.:0.04713 1st Qu.:4.400 1st Qu.:27.93 1st Qu.: 8.935
## Median :0.07840 Median :4.620 Median :30.67 Median :10.666
## Mean :0.08332 Mean :4.594 Mean :30.70 Mean :10.600
## 3rd Qu.:0.11800 3rd Qu.:4.870 3rd Qu.:33.75 3rd Qu.:12.297
## Max. :0.17460 Max. :5.180 Max. :39.05 Max. :16.475
## wbf.aud. stroke.amplitude wing.velocity totalMass
## Min. :163.9 Min. : 96.88 Min. :1.208 Min. :0.1062
## 1st Qu.:173.0 1st Qu.:114.16 1st Qu.:1.351 1st Qu.:0.1808
## Median :179.1 Median :125.96 Median :1.427 Median :0.2191
## Mean :178.9 Mean :123.31 Mean :1.465 Mean :0.2245
## 3rd Qu.:185.3 3rd Qu.:133.61 3rd Qu.:1.588 3rd Qu.:0.2646
## Max. :194.7 Max. :144.78 Max. :1.843 Max. :0.3550
## percLoad wingVel2 strokeAmp2
## Min. :1.143 Min. :7.978e-05 Min. : 0.0093
## 1st Qu.:1.321 1st Qu.:3.980e-03 1st Qu.: 42.9577
## Median :1.576 Median :1.444e-02 Median : 99.2962
## Mean :1.600 Mean :2.509e-02 Mean :148.2959
## 3rd Qu.:1.863 3rd Qu.:3.750e-02 3rd Qu.:188.1689
## Max. :2.180 Max. :1.423e-01 Max. :698.3210
## wbf2 totalMass2 percLoad2
## Min. : 0.00127 Min. :0.01128 Min. :0.0003255
## 1st Qu.: 6.59393 1st Qu.:0.03270 1st Qu.:0.0353246
## Median : 38.46046 Median :0.04801 Median :0.0764175
## Mean : 57.11835 Mean :0.05374 Mean :0.0964447
## 3rd Qu.: 77.70129 3rd Qu.:0.07003 3rd Qu.:0.1477398
## Max. :250.38161 Max. :0.12602 Max. :0.3363085
## wingVel_cent strokeAmp_cent wbf_cent
## Min. :-0.25706 Min. :-26.426 Min. :-15.074
## 1st Qu.:-0.11484 1st Qu.: -9.145 1st Qu.: -5.972
## Median :-0.03848 Median : 2.656 Median : 0.185
## Mean : 0.00000 Mean : 0.000 Mean : 0.000
## 3rd Qu.: 0.12225 3rd Qu.: 10.307 3rd Qu.: 6.394
## Max. : 0.37723 Max. : 21.473 Max. : 15.823
## totalMass_cent percLoad_cent
## Min. :-0.118253 Min. :-0.45679
## 1st Qu.:-0.043628 1st Qu.:-0.27912
## Median :-0.005353 Median :-0.02372
## Mean : 0.000000 Mean : 0.00000
## 3rd Qu.: 0.040172 3rd Qu.: 0.26297
## Max. : 0.130547 Max. : 0.57992
ggplot(bdta, aes(x = percLoad, fill = Treatment)) +
geom_histogram()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# seems like percent load and trt are very related
# visualize % load vs total mass
ggplot(bdta, aes(x = percLoad, y = totalMass, color = Treatment)) +
geom_point()
# seem to be strongly associated
# see number of observations per bee
table(bdta$Bee.ID) # each bee has two observations
##
## E01 E03 E05 E09 E10 E11 E12 E16 E20 E24 E26 E28 E30 E31 E32 E33 E35 E36
## 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## E37 E39 E41 E42 E44 E45 E49 E50 E53 E54 E56 E59
## 2 2 2 2 2 2 2 2 2 2 2 2
# get number of unque bees
length(table(bdta$Bee.ID))
## [1] 30
# get treatment orders
trtOrders_loadedSecond <- sapply(X = unique(bdta$Bee.ID), FUN = function(x){
tmp = bdta[bdta$Bee.ID == x, c("Treatment.order", "Treatment")]
if("2_L" %in% paste(tmp[,1], tmp[,2], sep = "_")){
loadedSecond = TRUE
}
else loadedSecond = FALSE
return(loadedSecond)
})
table(trtOrders_loadedSecond)
## trtOrders_loadedSecond
## FALSE
## 30
# visualize scatterplot
car::scatterplotMatrix(bdta[, 4:13])
# VIF is high among the bee size predictors
# note: this function is almost the same as car::vif
vif.mer(lmer(av.resp..CO2.mL.hr. ~ Treatment + Treatment.order + Mstarved + IT.Span + totalMass + single.wing.area + percLoad + (1|Bee.ID), data = bdta))
## TreatmentL Treatment.order2 Mstarved IT.Span
## 16.591645 1.270430 23.112962 6.658432
## totalMass single.wing.area percLoad
## 32.931859 10.903672 38.415618
car::scatterplotMatrix(bdta[, c("Mstarved", "IT.Span", "single.wing.area")])
# principle components
aa = prcomp(bdta[, c("Mstarved", "IT.Span", "single.wing.area")], center = TRUE, scale = TRUE)
summary(aa) # 1st pc explains ~95% of the variance in the three measurements of size
## Importance of components:
## PC1 PC2 PC3
## Standard deviation 1.6852 0.33843 0.21323
## Proportion of Variance 0.9467 0.03818 0.01516
## Cumulative Proportion 0.9467 0.98484 1.00000
biplot(aa) # shows that all three size measurement are correlated
# note, I changed the signs of the predictions so that higher PC1 values
# correspond to bigger bees
p1 = -predict(aa)[,1]
# add PC1 scores to dataset
bdta$size_pc1 = p1
# show scatterplot matrix to see correlations among size predictors
car::scatterplotMatrix(bdta[, c("Mstarved", "IT.Span", "single.wing.area", "size_pc1", "percLoad")])
# check VIF one more time
# VIF is high among the bee size predictors
vif.mer(lmer(av.resp..CO2.mL.hr. ~ Treatment + Treatment.order + size_pc1 + percLoad + totalMass + (1|Bee.ID), data = bdta))
## TreatmentL Treatment.order2 size_pc1 percLoad
## 16.643643 1.258990 7.448179 30.631628
## totalMass
## 26.790114
# remove treatment
vif.mer(lmer(av.resp..CO2.mL.hr. ~ Treatment.order + size_pc1 + percLoad + totalMass + (1|Bee.ID), data = bdta))
## Treatment.order2 size_pc1 percLoad totalMass
## 1.014834 7.698050 20.910198 25.813622
# looks like we cannot have total mass and percentLoad in the same model, according to VIF
# remove total mass
vif.mer(lmer(av.resp..CO2.mL.hr. ~ Treatment.order + size_pc1 + percLoad + (1|Bee.ID), data = bdta))
## Treatment.order2 size_pc1 percLoad
## 1.003811 1.004814 1.008624
# make a full model with all two-way interactions
m1 <- lmer(av.resp..CO2.mL.hr. ~ (size_pc1 + Treatment.order + percLoad_cent)^2 +
percLoad2 +
+ (1|Bee.ID), data = bdta)
summary(m1)
## Linear mixed model fit by REML ['lmerMod']
## Formula:
## av.resp..CO2.mL.hr. ~ (size_pc1 + Treatment.order + percLoad_cent)^2 +
## percLoad2 + +(1 | Bee.ID)
## Data: bdta
##
## REML criterion at convergence: 182.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.7335 -0.4895 0.1150 0.4689 2.3779
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 1.4520 1.2050
## Residual 0.5021 0.7086
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 11.23228 0.27961 40.17
## size_pc1 1.28130 0.15748 8.14
## Treatment.order2 -0.28829 0.18935 -1.52
## percLoad_cent 4.90035 0.70897 6.91
## percLoad2 -5.10731 1.38346 -3.69
## size_pc1:Treatment.order2 -0.17569 0.12108 -1.45
## size_pc1:percLoad_cent 0.09533 0.21461 0.44
## Treatment.order2:percLoad_cent -1.19963 1.15631 -1.04
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2 prcLd_ prcLd2 s_1:T. s_1:L_
## size_pc1 -0.001
## Trtmnt.rdr2 -0.253 0.032
## percLod_cnt -0.042 0.233 0.066
## percLoad2 -0.392 -0.049 -0.204 -0.043
## sz_pc1:Tr.2 -0.049 -0.410 -0.089 -0.325 0.219
## sz_pc1:prL_ -0.030 -0.105 -0.205 -0.146 0.331 0.289
## Trtmnt.2:L_ 0.071 -0.207 -0.029 -0.886 -0.035 0.284 0.159
### LRT's for interactions
m2.0 <- update(m1, .~. - size_pc1:percLoad_cent)
anova(m1, m2.0) # likelihood ratio test for interaction of treatment order
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m2.0: av.resp..CO2.mL.hr. ~ size_pc1 + Treatment.order + percLoad_cent +
## m2.0: percLoad2 + (1 | Bee.ID) + size_pc1:Treatment.order + Treatment.order:percLoad_cent
## m1: av.resp..CO2.mL.hr. ~ (size_pc1 + Treatment.order + percLoad_cent)^2 +
## m1: percLoad2 + +(1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2.0 9 195.92 214.77 -88.961 177.92
## m1 10 197.68 218.62 -88.841 177.68 0.2397 1 0.6244
summary(m2.0)
## Linear mixed model fit by REML ['lmerMod']
## Formula:
## av.resp..CO2.mL.hr. ~ size_pc1 + Treatment.order + percLoad_cent +
## percLoad2 + (1 | Bee.ID) + size_pc1:Treatment.order + Treatment.order:percLoad_cent
## Data: bdta
##
## REML criterion at convergence: 181.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.72009 -0.52773 0.07283 0.48826 2.35677
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 1.4603 1.2084
## Residual 0.4862 0.6973
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 11.2371 0.2783 40.37
## size_pc1 1.2881 0.1562 8.24
## Treatment.order2 -0.2709 0.1824 -1.49
## percLoad_cent 4.9339 0.6939 7.11
## percLoad2 -5.3211 1.2854 -4.14
## size_pc1:Treatment.order2 -0.1907 0.1141 -1.67
## Treatment.order2:percLoad_cent -1.2596 1.1316 -1.11
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2 prcLd_ prcLd2 s_1:T.
## size_pc1 -0.004
## Trtmnt.rdr2 -0.262 0.010
## percLod_cnt -0.048 0.219 0.036
## percLoad2 -0.401 -0.014 -0.147 0.008
## sz_pc1:Tr.2 -0.041 -0.393 -0.031 -0.299 0.136
## Trtmnt.2:L_ 0.078 -0.192 0.004 -0.885 -0.096 0.253
m2.1 <- update(m2.0,.~. - Treatment.order:percLoad_cent)
anova(m2.0, m2.1)
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m2.1: av.resp..CO2.mL.hr. ~ size_pc1 + Treatment.order + percLoad_cent +
## m2.1: percLoad2 + (1 | Bee.ID) + size_pc1:Treatment.order
## m2.0: av.resp..CO2.mL.hr. ~ size_pc1 + Treatment.order + percLoad_cent +
## m2.0: percLoad2 + (1 | Bee.ID) + size_pc1:Treatment.order + Treatment.order:percLoad_cent
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2.1 8 195.15 211.90 -89.573 179.15
## m2.0 9 195.92 214.77 -88.961 177.92 1.2255 1 0.2683
summary(m2.1)
## Linear mixed model fit by REML ['lmerMod']
## Formula:
## av.resp..CO2.mL.hr. ~ size_pc1 + Treatment.order + percLoad_cent +
## percLoad2 + (1 | Bee.ID) + size_pc1:Treatment.order
## Data: bdta
##
## REML criterion at convergence: 185
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.69131 -0.46047 0.00547 0.49254 2.34469
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 1.5578 1.2481
## Residual 0.4642 0.6813
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 11.2638 0.2811 40.07
## size_pc1 1.2546 0.1562 8.03
## Treatment.order2 -0.2696 0.1782 -1.51
## percLoad_cent 4.2487 0.3162 13.44
## percLoad2 -5.4887 1.2518 -4.38
## size_pc1:Treatment.order2 -0.1589 0.1079 -1.47
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2 prcLd_ prcLd2
## size_pc1 0.010
## Trtmnt.rdr2 -0.254 0.011
## percLod_cnt 0.044 0.102 0.086
## percLoad2 -0.383 -0.032 -0.147 -0.165
## sz_pc1:Tr.2 -0.061 -0.348 -0.033 -0.167 0.167
m2.2 <- update(m2.1, .~. - size_pc1:Treatment.order)
anova(m2.1, m2.2) ## drop all interactions
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m2.2: av.resp..CO2.mL.hr. ~ size_pc1 + Treatment.order + percLoad_cent +
## m2.2: percLoad2 + (1 | Bee.ID)
## m2.1: av.resp..CO2.mL.hr. ~ size_pc1 + Treatment.order + percLoad_cent +
## m2.1: percLoad2 + (1 | Bee.ID) + size_pc1:Treatment.order
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2.2 7 195.54 210.21 -90.773 181.54
## m2.1 8 195.15 211.90 -89.573 179.15 2.3985 1 0.1215
summary(m2.2)
## Linear mixed model fit by REML ['lmerMod']
## Formula:
## av.resp..CO2.mL.hr. ~ size_pc1 + Treatment.order + percLoad_cent +
## percLoad2 + (1 | Bee.ID)
## Data: bdta
##
## REML criterion at convergence: 184.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.6892 -0.4659 0.0722 0.3830 2.1677
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 1.5420 1.2418
## Residual 0.4856 0.6968
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 11.2373 0.2818 39.88
## size_pc1 1.1746 0.1462 8.03
## Treatment.order2 -0.2786 0.1822 -1.53
## percLoad_cent 4.1717 0.3187 13.09
## percLoad2 -5.1672 1.2616 -4.10
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2 prcLd_
## size_pc1 -0.012
## Trtmnt.rdr2 -0.261 -0.001
## percLod_cnt 0.035 0.049 0.081
## percLoad2 -0.385 0.029 -0.144 -0.141
# renaming model to simplify later typing
m2 <- m2.2
##### LRTs for main effects
## Treatment Order
m3 <- update(m2, .~. - Treatment.order)
anova(m2, m3, test = "Chi") # drop trt order (different than model with treatment)
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m3: av.resp..CO2.mL.hr. ~ size_pc1 + percLoad_cent + percLoad2 +
## m3: (1 | Bee.ID)
## m2: av.resp..CO2.mL.hr. ~ size_pc1 + Treatment.order + percLoad_cent +
## m2: percLoad2 + (1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m3 6 196.04 208.60 -92.019 184.04
## m2 7 195.54 210.21 -90.773 181.54 2.4925 1 0.1144
summary(m3)
## Linear mixed model fit by REML ['lmerMod']
## Formula: av.resp..CO2.mL.hr. ~ size_pc1 + percLoad_cent + percLoad2 +
## (1 | Bee.ID)
## Data: bdta
##
## REML criterion at convergence: 185.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.61825 -0.42513 -0.04217 0.52054 2.29988
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 1.5364 1.2395
## Residual 0.5073 0.7123
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 11.1235 0.2735 40.67
## size_pc1 1.1745 0.1465 8.02
## percLoad_cent 4.2122 0.3245 12.98
## percLoad2 -5.4323 1.2754 -4.26
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 prcLd_
## size_pc1 -0.013
## percLod_cnt 0.059 0.050
## percLoad2 -0.450 0.029 -0.131
# LRT for size
m4 <- update(m3, .~. - percLoad2)
anova(m3, m4) # keep percLoad2
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m4: av.resp..CO2.mL.hr. ~ size_pc1 + percLoad_cent + (1 | Bee.ID)
## m3: av.resp..CO2.mL.hr. ~ size_pc1 + percLoad_cent + percLoad2 +
## m3: (1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m4 5 208.25 218.73 -99.127 198.25
## m3 6 196.04 208.60 -92.019 184.04 14.216 1 0.0001629 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(m3)
## Linear mixed model fit by REML ['lmerMod']
## Formula: av.resp..CO2.mL.hr. ~ size_pc1 + percLoad_cent + percLoad2 +
## (1 | Bee.ID)
## Data: bdta
##
## REML criterion at convergence: 185.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.61825 -0.42513 -0.04217 0.52054 2.29988
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 1.5364 1.2395
## Residual 0.5073 0.7123
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 11.1235 0.2735 40.67
## size_pc1 1.1745 0.1465 8.02
## percLoad_cent 4.2122 0.3245 12.98
## percLoad2 -5.4323 1.2754 -4.26
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 prcLd_
## size_pc1 -0.013
## percLod_cnt 0.059 0.050
## percLoad2 -0.450 0.029 -0.131
# LRT for Treatment (load)
m5 <- update(m3, .~. - size_pc1)
anova(m3, m5) # keep size
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m5: av.resp..CO2.mL.hr. ~ percLoad_cent + percLoad2 + (1 | Bee.ID)
## m3: av.resp..CO2.mL.hr. ~ size_pc1 + percLoad_cent + percLoad2 +
## m3: (1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m5 5 229.84 240.31 -109.920 219.84
## m3 6 196.04 208.60 -92.019 184.04 35.801 1 2.185e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m6 <- update(m3, .~. - percLoad_cent)
anova(m3, m6) # keep squared term
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m6: av.resp..CO2.mL.hr. ~ size_pc1 + percLoad2 + (1 | Bee.ID)
## m3: av.resp..CO2.mL.hr. ~ size_pc1 + percLoad_cent + percLoad2 +
## m3: (1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m6 5 257.02 267.49 -123.510 247.02
## m3 6 196.04 208.60 -92.019 184.04 62.982 1 2.086e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# summarize final model for paper
summary(m3)
## Linear mixed model fit by REML ['lmerMod']
## Formula: av.resp..CO2.mL.hr. ~ size_pc1 + percLoad_cent + percLoad2 +
## (1 | Bee.ID)
## Data: bdta
##
## REML criterion at convergence: 185.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.61825 -0.42513 -0.04217 0.52054 2.29988
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 1.5364 1.2395
## Residual 0.5073 0.7123
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 11.1235 0.2735 40.67
## size_pc1 1.1745 0.1465 8.02
## percLoad_cent 4.2122 0.3245 12.98
## percLoad2 -5.4323 1.2754 -4.26
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 prcLd_
## size_pc1 -0.013
## percLod_cnt 0.059 0.050
## percLoad2 -0.450 0.029 -0.131
# write output
summary(m3)$coefficients
## Estimate Std. Error t value
## (Intercept) 11.123542 0.2734947 40.671875
## size_pc1 1.174509 0.1464512 8.019799
## percLoad_cent 4.212203 0.3245350 12.979195
## percLoad2 -5.432315 1.2754342 -4.259189
write.csv( summary(m3)$coefficients, file = "RespCoefs_percLoad.csv" )
# qq plot
qqnorm(resid(m3), main = "")
qqline(resid(m3)) # good
# residual plot
plot(fitted(m3), resid(m3), xlab = "fitted", ylab = "residuals")
abline(0,0)
# check cook's distance -- should be less than 1, so we're good
infl <- influence(m3, obs = TRUE)
plot(infl, which = 'cook')
# QQPlot for group-level effects
qqnorm(ranef(m3)$Bee.ID[[1]], main="Normal Q-Q plot for random effects")
qqline(ranef(m3)$Bee.ID[[1]]) # looks good
newdf <- data.frame(size_pc1 = 0, Treatment.order = factor(1), percLoad_cent = seq(min(bdta$percLoad_cent), max(bdta$percLoad_cent), length = 100))
newdf$percLoad2 <- newdf$percLoad_cent^2
preds1 <- predict(m3, newdf, re.form = NA)
# plot of prediction for average sized bee (line)
# raw data plotted as points
ggplot(bdta, aes(x = percLoad_cent + mean(bdta$percLoad) , y = av.resp..CO2.mL.hr.)) +
geom_point(alpha = 0.3) +
geom_line(data = newdf, aes(x = percLoad_cent + mean(bdta$percLoad), y = preds1)) +
labs(x = "Mass during flight (% of starved mass in g)")
# make a full model with all two-way interactions
m1 <- lmer(wbf.aud. ~ (size_pc1 + Treatment.order + percLoad_cent)^2 +
percLoad2 +
(1|Bee.ID), data = bdta)
summary(m1)
## Linear mixed model fit by REML ['lmerMod']
## Formula:
## wbf.aud. ~ (size_pc1 + Treatment.order + percLoad_cent)^2 + percLoad2 +
## (1 | Bee.ID)
## Data: bdta
##
## REML criterion at convergence: 326.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8304 -0.7064 0.1416 0.5616 1.5135
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 17.295 4.159
## Residual 9.876 3.143
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 183.2406 1.0745 170.53
## size_pc1 -2.9657 0.5916 -5.01
## Treatment.order2 -2.6627 0.8394 -3.17
## percLoad_cent 7.0660 2.9186 2.42
## percLoad2 -31.4157 6.0747 -5.17
## size_pc1:Treatment.order2 0.2089 0.5339 0.39
## size_pc1:percLoad_cent -0.5715 0.9461 -0.60
## Treatment.order2:percLoad_cent -0.7048 4.6314 -0.15
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2 prcLd_ prcLd2 s_1:T. s_1:L_
## size_pc1 0.005
## Trtmnt.rdr2 -0.294 0.038
## percLod_cnt -0.026 0.256 0.072
## percLoad2 -0.446 -0.066 -0.202 -0.081
## sz_pc1:Tr.2 -0.062 -0.476 -0.089 -0.311 0.228
## sz_pc1:prL_ -0.035 -0.121 -0.204 -0.146 0.332 0.286
## Trtmnt.2:L_ 0.061 -0.224 -0.034 -0.868 0.001 0.267 0.161
### LRT's for interactions
m2.0 <- update(m1, .~. - Treatment.order:percLoad_cent)
anova(m1, m2.0)
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m2.0: wbf.aud. ~ size_pc1 + Treatment.order + percLoad_cent + percLoad2 +
## m2.0: (1 | Bee.ID) + size_pc1:Treatment.order + size_pc1:percLoad_cent
## m1: wbf.aud. ~ (size_pc1 + Treatment.order + percLoad_cent)^2 + percLoad2 +
## m1: (1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2.0 9 362.18 381.03 -172.09 344.18
## m1 10 364.15 385.10 -172.08 344.15 0.0311 1 0.86
summary(m2.0)
## Linear mixed model fit by REML ['lmerMod']
## Formula:
## wbf.aud. ~ size_pc1 + Treatment.order + percLoad_cent + percLoad2 +
## (1 | Bee.ID) + size_pc1:Treatment.order + size_pc1:percLoad_cent
## Data: bdta
##
## REML criterion at convergence: 331.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8513 -0.6994 0.1341 0.5629 1.5305
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 16.921 4.113
## Residual 9.716 3.117
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 183.2499 1.0623 172.50
## size_pc1 -2.9860 0.5709 -5.23
## Treatment.order2 -2.6673 0.8321 -3.21
## percLoad_cent 6.6792 1.4349 4.65
## percLoad2 -31.4059 6.0246 -5.21
## size_pc1:Treatment.order2 0.2308 0.5103 0.45
## size_pc1:percLoad_cent -0.5481 0.9261 -0.59
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2 prcLd_ prcLd2 s_1:T.
## size_pc1 0.020
## Trtmnt.rdr2 -0.293 0.031
## percLod_cnt 0.054 0.128 0.086
## percLoad2 -0.447 -0.068 -0.202 -0.162
## sz_pc1:Tr.2 -0.082 -0.444 -0.082 -0.164 0.237
## sz_pc1:prL_ -0.045 -0.088 -0.201 -0.012 0.336 0.255
m2.1 <- update(m2.0,.~. - size_pc1:Treatment.order)
anova(m2.0, m2.1)
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m2.1: wbf.aud. ~ size_pc1 + Treatment.order + percLoad_cent + percLoad2 +
## m2.1: (1 | Bee.ID) + size_pc1:percLoad_cent
## m2.0: wbf.aud. ~ size_pc1 + Treatment.order + percLoad_cent + percLoad2 +
## m2.0: (1 | Bee.ID) + size_pc1:Treatment.order + size_pc1:percLoad_cent
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2.1 8 360.42 377.17 -172.21 344.42
## m2.0 9 362.18 381.03 -172.09 344.18 0.2349 1 0.6279
summary(m2.1)
## Linear mixed model fit by REML ['lmerMod']
## Formula:
## wbf.aud. ~ size_pc1 + Treatment.order + percLoad_cent + percLoad2 +
## (1 | Bee.ID) + size_pc1:percLoad_cent
## Data: bdta
##
## REML criterion at convergence: 332.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8110 -0.6782 0.1540 0.5820 1.4972
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 17.141 4.140
## Residual 9.392 3.065
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 183.2949 1.0536 173.96
## size_pc1 -2.8715 0.5122 -5.61
## Treatment.order2 -2.6345 0.8154 -3.23
## percLoad_cent 6.7941 1.3926 4.88
## percLoad2 -32.1204 5.7604 -5.58
## size_pc1:percLoad_cent -0.6559 0.8809 -0.74
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2 prcLd_ prcLd2
## size_pc1 -0.018
## Trtmnt.rdr2 -0.298 -0.006
## percLod_cnt 0.041 0.062 0.074
## percLoad2 -0.437 0.042 -0.189 -0.128
## sz_pc1:prL_ -0.025 0.028 -0.187 0.031 0.294
m2.2 <- update(m2.1, .~. - size_pc1:percLoad_cent)
anova(m2.1, m2.2) ## drop all interactions
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m2.2: wbf.aud. ~ size_pc1 + Treatment.order + percLoad_cent + percLoad2 +
## m2.2: (1 | Bee.ID)
## m2.1: wbf.aud. ~ size_pc1 + Treatment.order + percLoad_cent + percLoad2 +
## m2.1: (1 | Bee.ID) + size_pc1:percLoad_cent
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2.2 7 359.05 373.71 -172.53 345.05
## m2.1 8 360.42 377.17 -172.21 344.42 0.6319 1 0.4267
summary(m2.2)
## Linear mixed model fit by REML ['lmerMod']
## Formula:
## wbf.aud. ~ size_pc1 + Treatment.order + percLoad_cent + percLoad2 +
## (1 | Bee.ID)
## Data: bdta
##
## REML criterion at convergence: 334.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8541 -0.7036 0.1281 0.5915 1.3890
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 17.163 4.143
## Residual 9.258 3.043
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 183.2771 1.0500 174.54
## size_pc1 -2.8608 0.5114 -5.59
## Treatment.order2 -2.7473 0.7953 -3.45
## percLoad_cent 6.8290 1.3823 4.94
## percLoad2 -30.8797 5.4673 -5.65
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2 prcLd_
## size_pc1 -0.017
## Trtmnt.rdr2 -0.307 -0.001
## percLod_cnt 0.041 0.061 0.081
## percLoad2 -0.448 0.035 -0.143 -0.144
# renaming model to simplify later typing
m2 <- m2.2
##### LRTs for main effects
## Treatment Order
m3 <- update(m2, .~. - Treatment.order)
anova(m2, m3, test = "Chi") # keep treatment order
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m3: wbf.aud. ~ size_pc1 + percLoad_cent + percLoad2 + (1 | Bee.ID)
## m2: wbf.aud. ~ size_pc1 + Treatment.order + percLoad_cent + percLoad2 +
## m2: (1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m3 6 368.10 380.67 -178.05 356.10
## m2 7 359.05 373.71 -172.53 345.05 11.05 1 0.0008867 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# LRT for size
m4 <- update(m2, .~. - size_pc1)
anova(m2, m4) # keep size_pc1
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m4: wbf.aud. ~ Treatment.order + percLoad_cent + percLoad2 + (1 |
## m4: Bee.ID)
## m2: wbf.aud. ~ size_pc1 + Treatment.order + percLoad_cent + percLoad2 +
## m2: (1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m4 6 379.62 392.19 -183.81 367.62
## m2 7 359.05 373.71 -172.53 345.05 22.571 1 2.025e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# LRT for Treatment (load)
m5 <- update(m2, .~. - percLoad_cent)
anova(m2, m5) # keep percLoad
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m5: wbf.aud. ~ size_pc1 + Treatment.order + percLoad2 + (1 | Bee.ID)
## m2: wbf.aud. ~ size_pc1 + Treatment.order + percLoad_cent + percLoad2 +
## m2: (1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m5 6 375.91 388.47 -181.95 363.91
## m2 7 359.05 373.71 -172.53 345.05 18.858 1 1.408e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m6 <- update(m2, .~. - percLoad2)
anova(m2, m6) # keep squared term
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m6: wbf.aud. ~ size_pc1 + Treatment.order + percLoad_cent + (1 |
## m6: Bee.ID)
## m2: wbf.aud. ~ size_pc1 + Treatment.order + percLoad_cent + percLoad2 +
## m2: (1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m6 6 378.98 391.54 -183.49 366.98
## m2 7 359.05 373.71 -172.53 345.05 21.927 1 2.832e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# summarize final model for paper
summary(m2)
## Linear mixed model fit by REML ['lmerMod']
## Formula:
## wbf.aud. ~ size_pc1 + Treatment.order + percLoad_cent + percLoad2 +
## (1 | Bee.ID)
## Data: bdta
##
## REML criterion at convergence: 334.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8541 -0.7036 0.1281 0.5915 1.3890
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 17.163 4.143
## Residual 9.258 3.043
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 183.2771 1.0500 174.54
## size_pc1 -2.8608 0.5114 -5.59
## Treatment.order2 -2.7473 0.7953 -3.45
## percLoad_cent 6.8290 1.3823 4.94
## percLoad2 -30.8797 5.4673 -5.65
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2 prcLd_
## size_pc1 -0.017
## Trtmnt.rdr2 -0.307 -0.001
## percLod_cnt 0.041 0.061 0.081
## percLoad2 -0.448 0.035 -0.143 -0.144
# write output
summary(m2)$coefficients
## Estimate Std. Error t value
## (Intercept) 183.277063 1.0500400 174.542933
## size_pc1 -2.860783 0.5114480 -5.593497
## Treatment.order2 -2.747317 0.7952937 -3.454469
## percLoad_cent 6.829045 1.3823068 4.940325
## percLoad2 -30.879690 5.4672717 -5.648099
write.csv( summary(m2)$coefficients, file = "FreqCoefs_percLoad.csv" )
# qq plot
qqnorm(resid(m2), main = "")
qqline(resid(m2)) # good
# check influence
plot(influence(m2, obs = TRUE ), which = 'cook')
# residual plot
plot(fitted(m2), residuals(m2, "deviance"), xlab = "fitted", ylab = "residuals")
abline(0,0)
# QQPlot for group-level effects
qqnorm(ranef(m2)$Bee.ID[[1]], main="Normal Q-Q plot for random effects")
qqline(ranef(m2)$Bee.ID[[1]]) # looks good
newdf <- data.frame(size_pc1 = 0, Treatment.order = factor(1), percLoad_cent = seq(min(bdta$percLoad_cent), max(bdta$percLoad_cent), length = 100))
newdf$percLoad2 <- newdf$percLoad_cent^2
preds1 <- predict(m2, newdf, re.form = NA) #predict for mean of random effect
# plot of prediction for average sized bee (line)
# raw data plotted as points
ggplot(bdta, aes(x = percLoad_cent + mean(bdta$percLoad) , y = wbf.aud.)) +
geom_point(alpha = 0.3) +
geom_line(data = newdf, aes(x = percLoad_cent + mean(bdta$percLoad), y = preds1)) +
labs(x = "Mass during flight (% of starved mass in g)")
# make a full model with all two-way interactions
m1 <- lmer(stroke.amplitude ~ (size_pc1 + Treatment.order + percLoad_cent)^2 +
percLoad2 +
(1|Bee.ID), data = bdta)
summary(m1)
## Linear mixed model fit by REML ['lmerMod']
## Formula:
## stroke.amplitude ~ (size_pc1 + Treatment.order + percLoad_cent)^2 +
## percLoad2 + (1 | Bee.ID)
## Data: bdta
##
## REML criterion at convergence: 352.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.09359 -0.52135 0.01418 0.41271 2.08441
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 13.01 3.606
## Residual 24.48 4.948
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 124.1538 1.3504 91.94
## size_pc1 1.9451 0.7042 2.76
## Treatment.order2 -0.4425 1.3197 -0.34
## percLoad_cent 39.0584 3.8975 10.02
## percLoad2 -9.2460 9.3003 -0.99
## size_pc1:Treatment.order2 -1.3443 0.8311 -1.62
## size_pc1:percLoad_cent -3.1048 1.4603 -2.13
## Treatment.order2:percLoad_cent -5.6584 5.7085 -0.99
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2 prcLd_ prcLd2 s_1:T. s_1:L_
## size_pc1 0.025
## Trtmnt.rdr2 -0.374 0.049
## percLod_cnt 0.026 0.292 0.087
## percLoad2 -0.540 -0.105 -0.198 -0.167
## sz_pc1:Tr.2 -0.088 -0.611 -0.087 -0.284 0.242
## sz_pc1:prL_ -0.040 -0.151 -0.199 -0.147 0.329 0.278
## Trtmnt.2:L_ 0.016 -0.244 -0.047 -0.821 0.086 0.237 0.172
### LRT's for interactions
m2.0 <- update(m1, .~. - Treatment.order:percLoad_cent )
anova(m1, m2.0)
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m2.0: stroke.amplitude ~ size_pc1 + Treatment.order + percLoad_cent +
## m2.0: percLoad2 + (1 | Bee.ID) + size_pc1:Treatment.order + size_pc1:percLoad_cent
## m1: stroke.amplitude ~ (size_pc1 + Treatment.order + percLoad_cent)^2 +
## m1: percLoad2 + (1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2.0 9 394.21 413.05 -188.10 376.21
## m1 10 395.22 416.16 -187.61 375.22 0.9852 1 0.3209
summary(m2.0)
## Linear mixed model fit by REML ['lmerMod']
## Formula: stroke.amplitude ~ size_pc1 + Treatment.order + percLoad_cent +
## percLoad2 + (1 | Bee.ID) + size_pc1:Treatment.order + size_pc1:percLoad_cent
## Data: bdta
##
## REML criterion at convergence: 358.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.06103 -0.56098 -0.01219 0.50169 2.15159
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 14.18 3.765
## Residual 23.65 4.864
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 124.1660 1.3486 92.07
## size_pc1 1.7778 0.6856 2.59
## Treatment.order2 -0.5058 1.2961 -0.39
## percLoad_cent 35.9655 2.1941 16.39
## percLoad2 -8.3350 9.1404 -0.91
## size_pc1:Treatment.order2 -1.1503 0.7940 -1.45
## size_pc1:percLoad_cent -2.8411 1.4174 -2.00
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2 prcLd_ prcLd2 s_1:T.
## size_pc1 0.029
## Trtmnt.rdr2 -0.367 0.038
## percLod_cnt 0.068 0.163 0.085
## percLoad2 -0.536 -0.086 -0.196 -0.168
## sz_pc1:Tr.2 -0.094 -0.575 -0.078 -0.162 0.229
## sz_pc1:prL_ -0.044 -0.112 -0.195 -0.011 0.322 0.248
m2.1 <- update(m2.0,.~. - size_pc1:Treatment.order)
anova(m2.0, m2.1)
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m2.1: stroke.amplitude ~ size_pc1 + Treatment.order + percLoad_cent +
## m2.1: percLoad2 + (1 | Bee.ID) + size_pc1:percLoad_cent
## m2.0: stroke.amplitude ~ size_pc1 + Treatment.order + percLoad_cent +
## m2.0: percLoad2 + (1 | Bee.ID) + size_pc1:Treatment.order + size_pc1:percLoad_cent
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2.1 8 394.61 411.36 -189.3 378.61
## m2.0 9 394.21 413.05 -188.1 376.21 2.4002 1 0.1213
summary(m2.1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: stroke.amplitude ~ size_pc1 + Treatment.order + percLoad_cent +
## percLoad2 + (1 | Bee.ID) + size_pc1:percLoad_cent
## Data: bdta
##
## REML criterion at convergence: 362.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.02623 -0.51313 -0.04168 0.51460 2.13594
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 13.6 3.688
## Residual 24.7 4.970
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 123.9899 1.3564 91.41
## size_pc1 1.2045 0.5603 2.15
## Treatment.order2 -0.6514 1.3201 -0.49
## percLoad_cent 35.3986 2.2086 16.03
## percLoad2 -5.3965 9.0698 -0.59
## size_pc1:percLoad_cent -2.3436 1.4009 -1.67
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2 prcLd_ prcLd2
## size_pc1 -0.031
## Trtmnt.rdr2 -0.382 -0.008
## percLod_cnt 0.055 0.089 0.073
## percLoad2 -0.536 0.059 -0.183 -0.137
## sz_pc1:prL_ -0.021 0.040 -0.181 0.031 0.280
m2.2 <- update(m2.1, .~. - size_pc1:percLoad_cent)
anova(m2.1, m2.2) ## drop all interactions
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m2.2: stroke.amplitude ~ size_pc1 + Treatment.order + percLoad_cent +
## m2.2: percLoad2 + (1 | Bee.ID)
## m2.1: stroke.amplitude ~ size_pc1 + Treatment.order + percLoad_cent +
## m2.1: percLoad2 + (1 | Bee.ID) + size_pc1:percLoad_cent
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2.2 7 395.66 410.32 -190.83 381.66
## m2.1 8 394.61 411.36 -189.30 378.61 3.0552 1 0.08048 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(m2.2)
## Linear mixed model fit by REML ['lmerMod']
## Formula: stroke.amplitude ~ size_pc1 + Treatment.order + percLoad_cent +
## percLoad2 + (1 | Bee.ID)
## Data: bdta
##
## REML criterion at convergence: 367.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.89123 -0.52057 -0.02436 0.57562 2.02535
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 13.86 3.723
## Residual 25.63 5.062
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 123.9447 1.3781 89.94
## size_pc1 1.2414 0.5677 2.19
## Treatment.order2 -1.0513 1.3225 -0.79
## percLoad_cent 35.5007 2.2478 15.79
## percLoad2 -1.1743 8.8665 -0.13
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2 prcLd_
## size_pc1 -0.031
## Trtmnt.rdr2 -0.393 -0.001
## percLod_cnt 0.056 0.088 0.080
## percLoad2 -0.553 0.051 -0.140 -0.152
# renaming model to simplify later typing
m2 <- m2.2
##### LRTs for main effects
## size
m3 <- update(m2, .~. - percLoad2)
anova(m2, m3, test = "Chi") # drop squared term
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m3: stroke.amplitude ~ size_pc1 + Treatment.order + percLoad_cent +
## m3: (1 | Bee.ID)
## m2: stroke.amplitude ~ size_pc1 + Treatment.order + percLoad_cent +
## m2: percLoad2 + (1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m3 6 393.68 406.24 -190.84 381.68
## m2 7 395.66 410.32 -190.83 381.66 0.0174 1 0.8949
summary(m3)
## Linear mixed model fit by REML ['lmerMod']
## Formula: stroke.amplitude ~ size_pc1 + Treatment.order + percLoad_cent +
## (1 | Bee.ID)
## Data: bdta
##
## REML criterion at convergence: 373.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.87059 -0.50457 -0.03627 0.58077 2.02395
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 14.31 3.782
## Residual 24.73 4.973
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 123.8429 1.1414 108.50
## size_pc1 1.2463 0.5668 2.20
## Treatment.order2 -1.0742 1.2864 -0.84
## percLoad_cent 35.5004 2.1858 16.24
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2
## size_pc1 -0.003
## Trtmnt.rdr2 -0.564 0.006
## percLod_cnt -0.034 0.096 0.060
# trt order
m4 <- update(m3, .~. - Treatment.order)
anova(m3, m4) # drop treatment order
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m4: stroke.amplitude ~ size_pc1 + percLoad_cent + (1 | Bee.ID)
## m3: stroke.amplitude ~ size_pc1 + Treatment.order + percLoad_cent +
## m3: (1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m4 5 392.42 402.89 -191.21 382.42
## m3 6 393.68 406.24 -190.84 381.68 0.7374 1 0.3905
summary(m4)
## Linear mixed model fit by REML ['lmerMod']
## Formula: stroke.amplitude ~ size_pc1 + percLoad_cent + (1 | Bee.ID)
## Data: bdta
##
## REML criterion at convergence: 376.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.96941 -0.48117 -0.05909 0.48300 1.94459
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 14.53 3.811
## Residual 24.42 4.941
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 123.3058 0.9440 130.62
## size_pc1 1.2495 0.5675 2.20
## percLoad_cent 35.6293 2.1693 16.42
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1
## size_pc1 0.000
## percLod_cnt 0.000 0.095
# LRT for size
m5 <- update(m4, .~. - size_pc1)
anova(m4, m5) # keep size
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m5: stroke.amplitude ~ percLoad_cent + (1 | Bee.ID)
## m4: stroke.amplitude ~ size_pc1 + percLoad_cent + (1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m5 4 395.24 403.62 -193.62 387.24
## m4 5 392.42 402.89 -191.21 382.42 4.828 1 0.028 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m6 <- update(m4, .~. - percLoad_cent)
anova(m4, m6) # keep percent load
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m6: stroke.amplitude ~ size_pc1 + (1 | Bee.ID)
## m4: stroke.amplitude ~ size_pc1 + percLoad_cent + (1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m6 4 478.07 486.45 -235.04 470.07
## m4 5 392.42 402.89 -191.21 382.42 87.657 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# summarize final model for paper
summary(m4)
## Linear mixed model fit by REML ['lmerMod']
## Formula: stroke.amplitude ~ size_pc1 + percLoad_cent + (1 | Bee.ID)
## Data: bdta
##
## REML criterion at convergence: 376.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.96941 -0.48117 -0.05909 0.48300 1.94459
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 14.53 3.811
## Residual 24.42 4.941
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 123.3058 0.9440 130.62
## size_pc1 1.2495 0.5675 2.20
## percLoad_cent 35.6293 2.1693 16.42
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1
## size_pc1 0.000
## percLod_cnt 0.000 0.095
# write output
summary(m4)$coefficients
## Estimate Std. Error t value
## (Intercept) 123.305797 0.9440279 130.616690
## size_pc1 1.249463 0.5674547 2.201872
## percLoad_cent 35.629329 2.1693080 16.424283
write.csv( summary(m4)$coefficients, file = "AmpCoefs_percLoad.csv" )
# qq plot
qqnorm(resid(m4), main = "")
qqline(resid(m4)) # ok
# residual plot
plot(fitted(m4), resid(m3), xlab = "fitted", ylab = "residuals")
abline(0,0)
# QQPlot for group-level effects
qqnorm(ranef(m4)$Bee.ID[[1]], main="Normal Q-Q plot for random effects")
qqline(ranef(m4)$Bee.ID[[1]]) # looks good
preds1 <- predict(m4, newdf, re.form = NA) #predict for mean of random effect
# plot of prediction for average sized bee (line)
# raw data plotted as points
ggplot(bdta, aes(x = percLoad_cent + mean(bdta$percLoad) , y = stroke.amplitude)) +
geom_point(alpha = 0.3) +
geom_line(data = newdf, aes(x = percLoad_cent + mean(bdta$percLoad), y = preds1)) +
labs(x = "Mass during flight (% of starved mass in g)")
# make a full model with all two-way interactions
m1 <- lmer(wing.velocity ~ (size_pc1 + Treatment.order + percLoad_cent)^2 +
percLoad2 +
(1|Bee.ID), data = bdta)
summary(m1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: wing.velocity ~ (size_pc1 + Treatment.order + percLoad_cent)^2 +
## percLoad2 + (1 | Bee.ID)
## Data: bdta
##
## REML criterion at convergence: -76.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.96215 -0.44762 -0.04005 0.50444 2.11529
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 0.003968 0.06299
## Residual 0.006025 0.07762
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.48344 0.02182 67.98
## size_pc1 -0.04704 0.01148 -4.10
## Treatment.order2 -0.01763 0.02071 -0.85
## percLoad_cent -0.41560 0.06273 -6.62
## percLoad2 -0.05676 0.14666 -0.39
## size_pc1:Treatment.order2 0.01481 0.01306 1.13
## size_pc1:percLoad_cent 0.04749 0.02300 2.06
## Treatment.order2:percLoad_cent 0.05562 0.09323 0.60
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2 prcLd_ prcLd2 s_1:T. s_1:L_
## size_pc1 0.022
## Trtmnt.rdr2 -0.362 0.048
## percLod_cnt 0.017 0.288 0.085
## percLoad2 -0.527 -0.099 -0.199 -0.153
## sz_pc1:Tr.2 -0.085 -0.591 -0.087 -0.287 0.240
## sz_pc1:prL_ -0.040 -0.147 -0.200 -0.147 0.330 0.279
## Trtmnt.2:L_ 0.025 -0.242 -0.045 -0.829 0.072 0.241 0.170
### LRT's for interactions
m2.0 <- update(m1, .~. - Treatment.order:percLoad_cent )
anova(m1, m2.0)
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m2.0: wing.velocity ~ size_pc1 + Treatment.order + percLoad_cent +
## m2.0: percLoad2 + (1 | Bee.ID) + size_pc1:Treatment.order + size_pc1:percLoad_cent
## m1: wing.velocity ~ (size_pc1 + Treatment.order + percLoad_cent)^2 +
## m1: percLoad2 + (1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2.0 9 -101.513 -82.664 59.756 -119.51
## m1 10 -99.871 -78.927 59.935 -119.87 0.3579 1 0.5497
summary(m2.0)
## Linear mixed model fit by REML ['lmerMod']
## Formula: wing.velocity ~ size_pc1 + Treatment.order + percLoad_cent +
## percLoad2 + (1 | Bee.ID) + size_pc1:Treatment.order + size_pc1:percLoad_cent
## Data: bdta
##
## REML criterion at convergence: -79.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.02111 -0.46259 -0.01217 0.49529 2.12092
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 0.004027 0.06346
## Residual 0.005882 0.07670
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.48322 0.02168 68.42
## size_pc1 -0.04539 0.01109 -4.09
## Treatment.order2 -0.01704 0.02044 -0.83
## percLoad_cent -0.38478 0.03470 -11.09
## percLoad2 -0.06438 0.14470 -0.44
## size_pc1:Treatment.order2 0.01292 0.01253 1.03
## size_pc1:percLoad_cent 0.04506 0.02241 2.01
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2 prcLd_ prcLd2 s_1:T.
## size_pc1 0.028
## Trtmnt.rdr2 -0.359 0.038
## percLod_cnt 0.067 0.160 0.085
## percLoad2 -0.528 -0.084 -0.197 -0.167
## sz_pc1:Tr.2 -0.093 -0.561 -0.079 -0.162 0.230
## sz_pc1:prL_ -0.045 -0.110 -0.196 -0.011 0.324 0.249
m2.1 <- update(m2.0,.~. - size_pc1:Treatment.order)
anova(m2.0, m2.1)
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m2.1: wing.velocity ~ size_pc1 + Treatment.order + percLoad_cent +
## m2.1: percLoad2 + (1 | Bee.ID) + size_pc1:percLoad_cent
## m2.0: wing.velocity ~ size_pc1 + Treatment.order + percLoad_cent +
## m2.0: percLoad2 + (1 | Bee.ID) + size_pc1:Treatment.order + size_pc1:percLoad_cent
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2.1 8 -102.28 -85.525 59.140 -118.28
## m2.0 9 -101.51 -82.664 59.756 -119.51 1.2331 1 0.2668
summary(m2.1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: wing.velocity ~ size_pc1 + Treatment.order + percLoad_cent +
## percLoad2 + (1 | Bee.ID) + size_pc1:percLoad_cent
## Data: bdta
##
## REML criterion at convergence: -84.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.01405 -0.48637 0.04786 0.57874 2.09154
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 0.004062 0.06374
## Residual 0.005871 0.07662
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.485332 0.021598 68.77
## size_pc1 -0.038972 0.009194 -4.24
## Treatment.order2 -0.015368 0.020360 -0.75
## percLoad_cent -0.379047 0.034213 -11.08
## percLoad2 -0.099134 0.140706 -0.70
## size_pc1:percLoad_cent 0.039276 0.021687 1.81
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2 prcLd_ prcLd2
## size_pc1 -0.029
## Trtmnt.rdr2 -0.368 -0.008
## percLod_cnt 0.052 0.084 0.073
## percLoad2 -0.522 0.056 -0.184 -0.135
## sz_pc1:prL_ -0.022 0.038 -0.182 0.031 0.283
m2.2 <- update(m2.1, .~. - size_pc1:percLoad_cent)
anova(m2.1, m2.2) ## drop all interactions, though this is very close
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m2.2: wing.velocity ~ size_pc1 + Treatment.order + percLoad_cent +
## m2.2: percLoad2 + (1 | Bee.ID)
## m2.1: wing.velocity ~ size_pc1 + Treatment.order + percLoad_cent +
## m2.1: percLoad2 + (1 | Bee.ID) + size_pc1:percLoad_cent
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2.2 7 -100.72 -86.055 57.358 -114.72
## m2.1 8 -102.28 -85.525 59.140 -118.28 3.5643 1 0.05903 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(m2.2)
## Linear mixed model fit by REML ['lmerMod']
## Formula: wing.velocity ~ size_pc1 + Treatment.order + percLoad_cent +
## percLoad2 + (1 | Bee.ID)
## Data: bdta
##
## REML criterion at convergence: -87.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.89330 -0.56516 0.07513 0.49055 1.93807
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 0.004108 0.06410
## Residual 0.006188 0.07866
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.486075 0.022031 67.45
## size_pc1 -0.039588 0.009322 -4.25
## Treatment.order2 -0.008667 0.020552 -0.42
## percLoad_cent -0.380720 0.035078 -10.85
## percLoad2 -0.169751 0.138432 -1.23
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2 prcLd_
## size_pc1 -0.029
## Trtmnt.rdr2 -0.381 -0.001
## percLod_cnt 0.054 0.084 0.080
## percLoad2 -0.540 0.048 -0.141 -0.150
# renaming model to simplify later typing
m2 <- m2.2
##### LRTs for main effects
## size
m3 <- update(m2, .~. - Treatment.order)
anova(m2, m3, test = "Chi") # drop trt order
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m3: wing.velocity ~ size_pc1 + percLoad_cent + percLoad2 + (1 | Bee.ID)
## m2: wing.velocity ~ size_pc1 + Treatment.order + percLoad_cent +
## m2: percLoad2 + (1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m3 6 -102.52 -89.954 57.260 -114.52
## m2 7 -100.72 -86.055 57.358 -114.72 0.1952 1 0.6587
summary(m3)
## Linear mixed model fit by REML ['lmerMod']
## Formula:
## wing.velocity ~ size_pc1 + percLoad_cent + percLoad2 + (1 | Bee.ID)
## Data: bdta
##
## REML criterion at convergence: -93.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.94992 -0.54108 0.04463 0.51763 1.87287
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 0.004216 0.06493
## Residual 0.005997 0.07744
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.482722 0.020255 73.20
## size_pc1 -0.039609 0.009329 -4.25
## percLoad_cent -0.379848 0.034460 -11.02
## percLoad2 -0.179915 0.135098 -1.33
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 prcLd_
## size_pc1 -0.031
## percLod_cnt 0.090 0.083
## percLoad2 -0.643 0.048 -0.140
# trt order
m4 <- update(m3, .~. - percLoad2)
anova(m3, m4) # drop percLoad
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m4: wing.velocity ~ size_pc1 + percLoad_cent + (1 | Bee.ID)
## m3: wing.velocity ~ size_pc1 + percLoad_cent + percLoad2 + (1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m4 5 -102.80 -92.331 56.401 -112.80
## m3 6 -102.52 -89.954 57.260 -114.52 1.7177 1 0.19
summary(m4)
## Linear mixed model fit by REML ['lmerMod']
## Formula: wing.velocity ~ size_pc1 + percLoad_cent + (1 | Bee.ID)
## Data: bdta
##
## REML criterion at convergence: -93.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.05997 -0.53298 -0.03019 0.59699 2.09586
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 0.003737 0.06113
## Residual 0.006384 0.07990
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.465370 0.015198 96.42
## size_pc1 -0.038985 0.009136 -4.27
## percLoad_cent -0.385121 0.035065 -10.98
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1
## size_pc1 0.000
## percLod_cnt 0.000 0.095
# LRT for size
m5 <- update(m4, .~. - size_pc1)
anova(m4, m5) # keep size
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m5: wing.velocity ~ percLoad_cent + (1 | Bee.ID)
## m4: wing.velocity ~ size_pc1 + percLoad_cent + (1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m5 4 -89.649 -81.271 48.824 -97.649
## m4 5 -102.803 -92.331 56.401 -112.803 15.154 1 9.908e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m6 <- update(m4, .~. - percLoad_cent)
anova(m4, m6) # keep percent load
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m6: wing.velocity ~ size_pc1 + (1 | Bee.ID)
## m4: wing.velocity ~ size_pc1 + percLoad_cent + (1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m6 4 -48.945 -40.567 28.472 -56.945
## m4 5 -102.803 -92.331 56.401 -112.803 55.858 1 7.791e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# summarize final model for paper
summary(m4)
## Linear mixed model fit by REML ['lmerMod']
## Formula: wing.velocity ~ size_pc1 + percLoad_cent + (1 | Bee.ID)
## Data: bdta
##
## REML criterion at convergence: -93.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.05997 -0.53298 -0.03019 0.59699 2.09586
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 0.003737 0.06113
## Residual 0.006384 0.07990
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.465370 0.015198 96.42
## size_pc1 -0.038985 0.009136 -4.27
## percLoad_cent -0.385121 0.035065 -10.98
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1
## size_pc1 0.000
## percLod_cnt 0.000 0.095
# write output
summary(m4)$coefficients
## Estimate Std. Error t value
## (Intercept) 1.46537020 0.015197520 96.421669
## size_pc1 -0.03898547 0.009135556 -4.267444
## percLoad_cent -0.38512098 0.035064635 -10.983174
write.csv( summary(m4)$coefficients, file = "WingVelocity_percLoad.csv" )
# qq plot
qqnorm(resid(m4), main = "")
qqline(resid(m4)) # good
# cook's distance
plot(influence(m4, obs = TRUE), which = 'cook')
# residual plot
plot(fitted(m4), resid(m4), xlab = "fitted", ylab = "residuals")
abline(0,0)
# QQPlot for group-level effects
qqnorm(ranef(m4)$Bee.ID[[1]], main="Normal Q-Q plot for random effects")
qqline(ranef(m4)$Bee.ID[[1]]) # looks good
preds1 <- predict(m4, newdf, re.form = NA) #predict for mean of random effect
# plot of prediction for average sized bee (line)
# raw data plotted as points
ggplot(bdta, aes(x = percLoad_cent + mean(bdta$percLoad) , y = wing.velocity)) +
geom_point(alpha = 0.3) +
geom_line(data = newdf, aes(x = percLoad_cent + mean(bdta$percLoad), y = preds1)) +
labs(x = "Mass during flight (% of starved mass in g)")
# qq plot
qqnorm(resid(mm6), main = "")
qqline(resid(mm6)) # good, but two outliers
# residual plot
plot(fitted(mm6), resid(mm6), xlab = "fitted", ylab = "residuals")
abline(0,0)
# QQPlot for group-level effects
qqnorm(ranef(mm6)$Bee.ID[[1]], main="Normal Q-Q plot for random effects")
qqline(ranef(mm6)$Bee.ID[[1]]) # looks good
# visualize model:
sjp.lmer(mm6, type = 'fe')
## Computing p-values via Kenward-Roger approximation. Use `p.kr = FALSE` if computation takes too long.
newdf <- data.frame(strokeAmp_cent_scaled = 0,
wbf_cent_scaled = 0,
percLoad_cent_scaled = seq(min(bdta$percLoad_cent_scaled),
max(bdta$percLoad_cent_scaled),
length.out = 100),
Treatment.order = factor(1),
size_pc1_scaled = 0,
strokeAmp2_scaled = 0,
wbf2_scaled = 0
)
newdf$percLoad2_scaled <- scale(newdf$percLoad_cent_scaled^2)
preds1 <- predict(mm6, newdf, re.form = NA)
# plot of prediction for average sized bee (line)
# raw data plotted as points
ggplot(bdta, aes(x = percLoad_cent_scaled , y = av.resp..CO2.mL.hr.)) +
geom_point(alpha = 0.3) +
geom_line(data = newdf, aes(x = percLoad_cent_scaled, y = preds1)) +
labs(x = "Mass during flight (% of starved mass in g)")
# visualize metabolic rate vs. wbf
newdf <- data.frame(strokeAmp_cent_scaled = 0,
wbf_cent_scaled = seq(min(bdta$wbf_cent_scaled),
max(bdta$wbf_cent_scaled),
length.out = 100),
percLoad_cent_scaled = 0,
percLoad2_scaled = 0,
Treatment.order = factor(1),
size_pc1_scaled = 0,
strokeAmp2_scaled = 0
)
newdf$wbf2_scaled <- scale(newdf$wbf_cent_scaled^2)
preds1 <- predict(mm6, newdf, re.form = NA)
# plot of prediction for average sized bee (line)
# raw data plotted as points
ggplot(bdta, aes(x = wbf_cent_scaled , y = av.resp..CO2.mL.hr.)) +
geom_point(alpha = 0.3) +
geom_line(data = newdf, aes(x = wbf_cent_scaled, y = preds1)) +
labs(x = "Wing beat frequency (centered and scaled)")
# visualize metabolic rate vs. stroke
newdf <- data.frame(
wbf_cent_scaled = 0,
wbf2_scaled = 0,
percLoad_cent_scaled = 0,
percLoad2_scaled = 0,
Treatment.order = factor(1),
size_pc1_scaled = 0,
strokeAmp_cent_scaled = seq(min(bdta$strokeAmp_cent_scaled),
max(bdta$strokeAmp_cent_scaled),
length.out = 100)
)
newdf$strokeAmp2_scaled <- scale(newdf$strokeAmp_cent_scaled^2)
preds1 <- predict(mm6, newdf, re.form = NA)
# plot of prediction for average sized bee (line)
# raw data plotted as points
ggplot(bdta, aes(x = strokeAmp_cent_scaled , y = av.resp..CO2.mL.hr.)) +
geom_point(alpha = 0.3) +
geom_line(data = newdf, aes(x = strokeAmp_cent_scaled, y = preds1)) +
labs(x = "Stroke Amplitude (centered and scaled)")
This model is fundamentally different from the previous one. It’s saying, when we account for differences in wingbeat freq, stroke amplitude, and bee size, we sill see an effect of percent load on metabolic rate.
OR, we could say that holding all other variables constant, an increase in wingbeat frequency is associated with an increase in metabolic rate.